Characterization of the material response in the granular ratcheting 
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The existence of a very special ratcheting regime has recently been reported in a granular packing 
subjected to cyclic loading [1]. In this state, the system accumulates a small permanent deformation 
after each cycle. After a short transient regime, the value of this permanent strain accumulation 
becomes independent on the number of cycles. We show that a characterization of the material 
response in this peculiar state is possible in terms of three simple macroscopic variables. They 
are defined that, they can be easily measured both in the experiments and in the simulations. We 
have carried out a thorough investigation of the micro- and macro-mechanical factors affecting these 
variables, by means of Molecular Dynamics simulations of a polydisperse disk packing, as a simple 
model system for granular material. Biaxial test boundary conditions with a periodically cycling 
load were implemented. The effect on the plastic response of the confining pressure, the deviatoric 
stress and the number of cycles has been investigated. The stiffness of the contacts and friction has 
been shown to play an important role in the overall response of the system. Specially elucidating is 
the influence of the particular hysteretical behavior in the stress-strain space on the accumulation 
of permanent strain and the energy dissipation. 
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I. INTRODUCTION 

Brownian motors, quantum ratchets or molecular 
pumps, all these machines operate under the same prin- 
ciple: The chaos of the micro-world cannot be avoided, 
but one can take advantage of it [2] . Nanoscale ratchet 
devices have been designed with the surprising property 
that they can extract work from the noise of thermal 
and quantum fluctuations [3]. Ratcheting is the mecha- 
nism behind molecular motors, which can use the chaotic 
Brownian motion to turn directionless energy into di- 
rected motion [4] . These lilliputian motors seem to be re- 
sponsible for many biological process, such as mechanical 
transport [5] or muscle contraction [6] . Apart from these 
fascinating machines, the ratchet effect has been used to 
describe economical or sociological processes where the 
intrinsic asymmetry in the system allow to rectify an un- 
biased input [7]. A ratchet-like effect is also the major 
cause of material deterioration due to cyclic stress load- 
ing, thermal or mechanical fluctuations [8-10]. Asymme- 
tries in foundations can produce tilting and eventual col- 
lapse of any structure due to ratcheting [11]. The tower 
of Pisa is a well documented case, where the tilt was ob- 
served from its construction in 1173 [12]. Pavement de- 
sign is another important field in which graded soils are 
used as supportive roadbed [8, 13-15]. The excitations 
that traffic imposes on the sub-layer produce deforma- 
tions in the granular material. These deformations are 
transmitted to the upper layers of the pavements, caus- 
ing its degradation or even its breakage. Cyclic load- 
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ing tests are extensively used in the investigation of the 
plastic response of unbound granular matter [13] . In 
these experiments, the material is subjected to a certain 
cyclic stress condition mimicking traffic. From a practi- 
cal point of view, the main question is whether the ma- 
terial accumulates plastic deformation in each cycle, or 
whether it adapts to the excitation reaching a shakedown 
state. Only materials in which the excitations shake down 
should be consequently used in pavement design. 

The use of simple models of granular materials allows 
the numerical solution of the dynamics. Discrete Ele- 
ment Methods (DEM) such as Molecular Dynamics (MD) 
[16-18] and Contact Dynamics (CD) [19, 20] have been 
in fact often successfully applied to the investigation of 
the elasto-plastic behavior of granular matter. Specially 
interesting from the physical point of view, is how the 
contact modelization affects the overall response [21, 22]. 
Recent MD results have shown the key role that sliding 
plays on the plastic deformation of a granular packing 
subjected to cyclic loading, and the existence of a range 
of values of the excitations for which a simple visco-elastic 
model of disks subjected to cyclic loading attains shake- 
down [1, 23]. Beyond the shakedown limit, two other pos- 
sible responses have also been identified; For very high 
loads, the material accumulates deformations at a rela- 
tively high constant rate, leading to an incremental col- 
lapse of the structure. For moderate loading intensities, 
the system undergoes an adaptation process in which 
the accumulation of deformation gradually decreases to 
a very low constant value. This post-compaction is asso- 
ciated to a relaxation of the dissipated energy per cycle, 
that progressively decreases to a constant value depen- 
dent on the imposed loading. In this final stage, there is 
a small but persistent accumulation of permanent strain, 
associated to a periodic behavior of the sliding contacts 
[1], which is called ratcheting regime. 
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Due to the non-lineality and the irreversibility of the 
behavior, cyclic loading is a rather complicated problem 
from the theoretical point of view. Elasto-plastic and hy- 
poplastic theories can account for the change in the in- 
cremental stiffness during loading and unloading phases, 
only if basic modifications are undertaken [24, 25]. In 
the case of elastoplasticity, the overall plastic behavior in 
the loading-unloading is obtained as the result of a com- 
bination of several yield surfaces [26] . In the hypoplastic 
theory the inter- granular strain is introduced to take into 
account the dependence of the response on the deforma- 
tion history [27]. Interestingly, a point of convergence 
of both theories has been established by the bounding 
surface elasto-plasticity [28]. This theory introduces a 
tiny elastic nucleus changing with the deformation, and 
describe the hysteresis by means of internal variables tak- 
ing into account the evolution of the microstructure. The 
characterization of such internal variables has been tra- 
ditionally done using structure tensors, measuring the 
fabric properties of the contact network [29]. There is 
numerical evidence that a single fabric tensor, measur- 
ing the anisotropy of the contact network can be used 
to characterize the resilient response [30]. But the de- 
scription of the plastic deformation requires to take into 
account the inherent decomposition of the contact net- 
work in sliding and non-sliding contacts [31]. The role of 
kinematical modes such as sliding and rolling has been 
also investigated to some extent for monotonic deforma- 
tion, but not for cyclic loading [16, 32, 33]. 

The final aim of this paper is the characterization of 
ratcheting response of a granular packing under cyclic 
loading. For this purpose, three macroscopic variables 
will be introduced. A simple DEM model will then be 
used to investigate the dependence of the material re- 
sponse on different macroscopic and microscopic vari- 
ables. From this investigation, we have found that our 
simple model is able to reproduce several behaviors ob- 
served in the experience, and microscopically justifies the 
use of popular empirical laws, like the k — 9 model. The 
main parameters of our model and the details of the MD 
simulations are presented in Section II. The ratcheting 
regime resulting in the biaxial test is described in Sect. 
III. In Sect. IV we decompose the strain response in its 
permanent and resilient components. We continue with 
an analysis of hysteresis in the plastic response, establish- 
ing in Sect. V, a direct relation between the particular 
shape of the stress-strain cycle and the dissipated energy 
per cycle. From this relationship it will be easy to ex- 
plain the observed dependence of the dissipated energy 
per cycle on the deviatoric stress. Results on the per- 
manent strain and the resilient parameters are presented 
for the different cases studied in Sees. VI & VII. The 
approach proposed here is basically empirical. The re- 
silient parameters will be therefore conveniently defined 
in terms of the recoverable deformation, as is usually 
done by experimentalists [34]. The dependence on the 
imposed stress is investigated, and the results are com- 
pared to predictions of resilient response models [35-38] . 



The influence of the friction and the stiffness at the con- 
tacts, main micro-mechanical parameters of the model, 
will also be determined. We finish in Section VIII with 
a discussion of the main conclusions of this work. 



II. MODEL 

In our visco-elastic 2D model, the grains are modeled 
by soft disks. The deformation that two grains suffer 
during the interaction is reproduced by letting the disks 
overlap. During the overlapping, a certain force f c is ex- 
erted at the contact point. This force can be decomposed 
in the following parts: 

r = r + r, w 

where f e and /" are the elastic and viscous contribution. 
The elastic part of the contact force is also decomposed 
as 

f e = /> c + f t H c . (2) 

The unit normal vector n c points in the direction of the 
vector connecting the center of mass of the two disks. 
The tangential vector t c is perpendicular to n c . The nor- 
mal elastic force is calculated as 

f e n = -k n A/L c , (3) 

where k n is the normal stiffness, A is the overlapping 
area and L c is a characteristic length of the contact. Our 
choice is L c = Ri + Rj . This normalization is necessary 
to be consistent in the units of force. 

The frictional force is calculated using an extension of 
the method proposed by Cundall-Strack [39]. An elastic 
force proportional to the elastic displacement is included 
at each contact 

ft = -k t Ax e t , (4) 

where k t is the tangential stiffness. The elastic displace- 
ment Axt is calculated as the time integral of the tan- 
gential velocity of the contact during the time where the 
elastic condition |/ t e | < /j,f% is satisfied. The sliding 
condition is imposed, keeping this force constant when 
|/ t e | = /j,f^. The straightforward calculation of this elas- 
tic displacement is given by the time integral starting at 
the beginning of the contact: 

Ax\ = f vt{t')Q{»t n -\m)dt', (5) 
Jo 

where is the Heaviside step function and v% denotes 
the tangential component of the relative velocity v c at 
the contact: 
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Vi 



Vj +uji x Ri 



a 3 x Rj. 



(6) 



Here Vi is the velocity and u5j is the angular velocity of 
the particles in contact. The branch vector Ri connects 
the center of mass of particle i to the point of application 
of the contact force. Replacing Eqs. (3) and (4) into (2) 
one obtains: 



P 



k t Ax e A c 



(7) 



Damping forces are included in order to allow rapid 
relaxation during the preparation of the sample, and to 
reduce the acoustic waves produced during the loading. 
These forces are calculated as 



P = -m( ln v c n n c + lt v c t t c ), 



(8) 



_1 the effective mass of the 
are the normal and tangen- 



bcing m = (l/m^ + 1/rrij 
disks in contact. n c and t 
tial unit vectors defined before, and 7„ and "f t are the 
coefficients of viscosity. These forces introduce time de- 
pendent effects during the loading. However, these effects 
can be arbitrarily reduced by increasing the loading time, 
as corresponds to the quasi-static approximation. 

The interaction of the disks with the walls is modeled 
by using a simple visco-elastic force: First, we allow the 
disks to penetrate the walls. Then we include a force 



P = - (M + lbm a v h ) 



n. 



(9) 



where 8 is the penetration length of the disk, n is the unit 
normal vector to the wall, and v b is the relative velocity 
of the disk with respect to the wall. 

The evolution of the position a?j and the orientation ipi 
of the particle i is governed by the equations of motion: 



c b 

im = Y,Ri*P + Y. % x ( 10 ) 

c b 

Here m, and 7j are the mass and moment of inertia 
of the disk. The first sum goes over all those particles 
in contact with this particle; the second one over all the 
forces given by the walls. The interparticle contact forces 
f c are given by replacing Eqs. (7) and (8) in Eq. (1). 

We use a fifth-order Gear predictor-corrector method 
for solving the equation of motion [40]. This algorithm 
consists of three steps. The first step predicts position 
and velocity of the particles by means of a Taylor expan- 
sion. The second step calculates the forces as a function 
of the predicted positions and velocities. The third step 
corrects the positions and velocities in order to optimize 



the stability of the algorithm. This method is much more 
efficient than the simple Euler approach or the Runge- 
Kutta method, especially for cyclic loading, where very 
high accuracy is required. 

The relevant contact parameters of this model are the 
normal stiffness at the contacts k n , the ratio of tangen- 
tial and normal stiffness kt/k n , the normal and tangential 
damping frequencies and the friction coefficient. In the 
quasi-static approximation, the results are independent 
of the frequency of the cyclic loading. The system is 
polydisperse, being the radii of the grains Gaussian dis- 
tributed with mean value of 1.0cm and variance of 0.36. 



III. ONSET OF GRANULAR RATCHETING 

In a biaxial experiment, the sample is subjected to a 
certain stress state characterized by the principal stresses 
o\ and ct 2 . In this case the stress space is therefore a 
plane, since the third component is zero, 03 = 0. In 
our simulations, the system is first homogeneously com- 
pressed with cTi = (72- After an equilibrium state under 
the pressure Po = 0-1 +°" 2 = has been reached, the 
vertical stress is quasi-statically changed: 



<r 2 (t) = Po 



cos 



2vrf 



(11) 



where t is the simulation time and to is the period of the 
loading. Note that Act, introduced in the last equation, 
is the maximum deviatoric stress measured in units of Po- 
In our approximation, it fully characterizes the intensity 
of the cyclic load imposed on the walls. 

Deformation appears in the sample due to the imposed 
excitations. The strain is the magnitude that charac- 
terizes the accumulation of permanent deformation in 
the sample. Among the different practical definitions 
of strain available [41], we have chosen here Cauchy's 
definition, which is basically the ratio of the new and 
the original length of the system. Let L be the original 
length of the sample in the principal direction i (i = x,y). 
The principal component of the strain tensor on this 
direction will then be: 



ei(t) = e u (t) = 



Li(t) — L Q 



(12) 



where Li is the length of the system in the principal di- 
rection i at the moment of the measurement. 

Different loading intensities will be exerted on the sam- 
ple by changing the value of Act. The reaction of the 
system will be characterized by the deviatoric permanent 
strain, 7, that is the difference between the strains in the 
principal directions: 



7 = e 2 - ei. 



(13) 
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FIG. 1: Typical stress-strain relation during cyclic loading. 
In the long-time behavior the response is given by a limit 
hysteresis loop. This is shown by the dashed line ( the loop 
here corresponds to N = 1000). In this simulation Act = 
0.14Po and Po = 10~ 4 fc„, where the normal contact stiffness 
is k„ = 2 ■ 10 6 N/m. The damping constants are defined in 
terms of the characteristic oscillation period t s = sjk n /p\ 2 
(in our case t s = 0.1414), where p is the density of the grains 
and A the mean radius of the disks composing the sample. 
The period of oscillation was taken long enough (to = 10 5 i s ), 
to be sure that we are in the quasi-static limit). 



The typical evolution of the permanent strain during 
the cyclic loading is shown in Fig. 1. The stress-strain 
relation consists of hysteresis loops. This hysteresis pro- 
duces an accumulation of deviatoric strain with the num- 
ber of cycles in addition to a progressive compaction, 
which is not shown there. After some decades of cycles, 
the accumulation of permanent deformation becomes lin- 
ear, as shown in Fig. 2. This strain rate remains constant 
for very large number of cycles, even when the volume ra- 
tio is very close to the saturation level. 

A micro-mechanical explanation of this linear accumu- 
lation of strain is provided by following the dynamics of 
the contact network. Although most of the contact forces 
of this network satisfy the elastic condition < jtt/ n , 
the strong heterogeneities produce a considerable amount 
of contacts reaching the sliding condition \f t \ = fif n dur- 
ing the compression. After a number of loading cycles, 
the contact network reaches a quasi-periodic behavior. 
In this regime, a fraction of the contacts reaches almost 
periodically the sliding condition, as shown in the inset of 
Fig. 2. In each load-unload transition there is an abrupt 
reduction of sliding contacts, which induces the typical 
discontinuity of the stiffness upon reversal of the load- 
ing. The load-unload asymmetry at each sliding contact 
makes it to slip the same amount and in the same di- 
rection during each loading cycle, leading to an overall 



ratcheting response. 

The contact behavior can be observed by embedding 
two points at each particle near to the contact area, and 
following their translation during each cycle. Their rela- 
tive displacements are calculated as s 1 = sq — s,.t>, where 
so is the displacement of the embedded point i and s rc i is 
the rigid body motion. This latter is given by the vector 
connecting the initial to the final position of the contact 
point. Note that s = when the two particles move as a 
rigid body. 

Figure 3 shows the displacement at the contacts during 
cycle N = 1000. Simulations show that, in this regime, 
this displacement field is almost constant after each cy- 
cle. There are two deformation modes resembling the 
mechanical ratchets: (i) At the sliding contacts the dis- 
placement vectors do not agree, so that there is a system- 
atic slip during each cycle which also leads to a constant 
frictional dissipation per cycle, (ii) At the non-sliding 
contacts the displacement vectors are almost the same 
for the two particles. 

Note from Fig. 3 that the distribution of this ratchets 
are not uniform, but they are localized in layers resem- 
bling shear bands. This kind of strain localization with 
intense rolling is typical in sheared granular materials 
[32, 42]. Fundamental differences are however observed 
between the cyclic loading response and the behavior un- 
der monotonic shear: The translation of each particle 
during the ratcheting regime is given by an almost con- 
stant displacement per cycle. On the other hand, the 
displacement of the particle during monotonic shear is 
rather chaotic, well described by an anomalous diffusion 
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FIG. 2: Cumulative permanent deformation against the num- 
ber of cycles (N). After the post-compaction regime the sys- 
tem accumulates permanent strain at a constant strain-rate. 
This is the so-called ratcheting regime, which emerges as a 
result of the periodicity of the sliding contacts. The inset 
precisely shows the fraction of the sliding contact versus time 
in this state. 
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FIG. 3: Displacement at the contacts during one cycle in the 
ratcheting. The arrows are proportional to the displacements 
s of the two material points at the contacts referred to the 
contact point. More details are found in the text. The figure 
is a snapshot of the simulation of Figure 1 for N = 1000. 



FIG. 4: Vortex formation as a consequence of the ratcheting of 
the particles. The arrows are proportional to the displacement 
of the particle after one cycle in the ratcheting regime. They 
are plotted at the center of the disks. The cycle is the same 
as the one shown in Figure 3. 



IV. MATERIAL RESPONSE TO CYCLIC 
LOADING 



[43]. 

Such systematic translation per cycle of the individ- 
ual grains in the ratcheting regime has a strong spatial 
correlation. This is shown in the displacement field of 
Fig. 4. The most salient feature here is the formation 
of vorticity cells, where a cluster of particles rotates as a 
whole. These vorticities survive during several hundred 
of cycles. This is contrary to the case of the simple shear, 
where the vorticities have a very short life-time [43]. It 
is interesting to see from Figures 3 and 4 the kinematic 
phase separation of the grains: (a) Grains organized in 
large vorticity cells, and (b) grains which accommodate 
the cells to make them more compatible with the imposed 
boundary conditions. Since such kinematical modes are 
linked with the a non-vanishing antisymmetric part of 
the displacement gradient, the strain tensor is not suffi- 
cient to provide a complete description of this convective 
motion during cyclic loading. An appropriate continuum 
description of ratcheting would require additional con- 
tinuum variables taking into account the vorticity and 
the gearing between the contacts. As in the case of the 
shear band formation, the Cosserat theory may be a good 
alternative [44] 



The existence of an elastic region in the deformation of 
granular materials implies that there is a finite region in 
the space of stress-states around the origin, in which the 
system reacts reversibly. Experiments and simulations 
show, however, that there is not such pure elastic behav- 
ior in a granular sample. Note, that this is not in con- 
tradiction with the existence of shakedown: A granular 
system may not accumulate any systematic permanent 
deformation after one loading cycle, but will always dis- 
sipate some energy because grain interactions are inher- 
ently inelastic. This is possible thanks to the additional 
energy supplied to the system by the external loading. 
In the particular case of our model, the system reaches 
a visco-elastic shakedown. In this limit state, the sys- 
tem dissipates some energy in each cycle and the overall 
behavior is not elastic, but the stress-strain cycle is still 
hysteretic (see Figure 5). Therefore we differentiate in 
cyclic loading between an elastic and a resilient deforma- 
tion of the sample. The latter implying that no perma- 
nent deformation has been accumulated after one cycle, 
while the first also implies the total absence of hysteresis 
or memory effects in the response. 

It has been recently shown that there is a broad range 
of values of Act for which a granular packing reacts to 
the imposed cyclic excitations by slowly deforming in a 
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ratcheting regime [1, 23]. This is a quasi-periodic state, 
macroscopically characterized by a constant strain rate 
and a conservation of the shape of the stress-strain cy- 
cle (see Figure 1). At the beginning of the loading pro- 
cess, the system suffers a re-arrangement of the sliding 
contacts, after which they start to behave periodically 
within the loading cycles. This post- compaction process 
is associated to a relaxation of the strain rate and also of 
the dissipated energy per cycle towards a constant value 
[23] . This stationary value of the strain rate fully deter- 
mines the macroscopic plastic response of the system in 
the ratcheting regime. At any stage of the experiment, 
the strain can therefore be decomposed in two well differ- 
entiated components. The irreversible plastic strain ac- 
cumulated after the end of the current cycle jp, and the 
recoverable resilient strain, -jr, accumulated along the 
cycle. In the ratcheting regime, the strain rate (A7/A./V) 
is approximately constant, while the latter deformation 
is well characterized by the resilient parameters: resilient 
modulus, Mr and the Poisson ratio £. The first param- 
eter, as it appears in Figure 5, is the ratio of the maxi- 
mum deviatoric stress and the corresponding deviatoric 
resilient strain: 



Ay/ AN 



Mu = 



Act 

1R ' 



(14) 



and quantifies the overall stiffness of the material. The 
Poisson ratio, correspondingly, is the ratio of the hori- 
zontal (ef ) and axial (ef ) resilient strains: 



ef 



C = ~e« 



(15) 



It is a measure for the isotropy of the deformation. The 
definition of ef and ef is similar to that in Eq. (12). 
They are both measured at the final stage of the load- 
ing, just before unloading starts. Similarly to Eq. (13), 
the resilient deviatoric strain is defined in terms of the 
resilient strains as y R — ef — ef . 

As a consequence of the quasi-static change of the 
stresses, all the relevant time dependence occurs in the 
system through the number of cycles N. Figure 6 shows 
the evolution of the resilient parameters from the simu- 
lations for different deviatoric stresses. For low excita- 
tions, the curves have already reached a plateau after a 
couple of cycles, implying that the values of ( and Mr do 
not apparently change as the number of cycles increases. 
In the initial post-compaction stage, the system accu- 
mulates more deviatoric strain in the horizontal direc- 
tion (perpendicular to the direction on which the cyclic 
load is applied), than it does in the final stage. This ex- 
plains why Poisson ratio decreases slightly in the first cy- 
cles. The resilient modulus increases, however, implying 
a higher stiffness of the system after the post-compaction. 
Although the dependence of the final values on the im- 
posed loading will be discussed in a latter section of this 
paper, it should now be remarked that the number of 
cycles needed for the system to reach a steady resilient 




FIG. 5: Sketch of the typical material reaction to cyclic load- 
ing in the granular ratcheting. After a post-compaction stage, 
the system accumulates permanent strain, 7p, at a constant 
strain-rate Ay /AN. The resilient modulus Mr is also indi- 
cated, as defined in Eq. (14). 



response increases as the imposed deviatoric stress is in- 
creased. This is clearly observed in case Act = 0.35 of 
the figure, where even after N = 1000 cycles, neither £ 
nor Mr have reached a stationary value. 

The peculiar behavior of the system in the ratchet- 
ing regime allows for the characterization of the defor- 
mation state of the system through the strain rate and 
the resilient parameters. It is therefore crucial to know 
the influence of the confining pressure and the deviatoric 
stress on these parameters. For a complete review on the 
macroscopic factors affecting the resilient response of a 
granular material and some of the models proposed to 
account for it, we recommend references [34] and [45]. 

To our knowledge, no systematic study has been car- 
ried out up to now elucidating the effect of the micro- 
scopic parameters of the system on the material reaction 
to cyclic loading, although they play an important role 
[46, 47]. Combe et al. have identified contact stiffness 
and friction as the relevant microscopic parameters in 
this limit. Inter-granular friction, in particular, appears 
then to be the dominating dissipative mechanism. The 
influence of contact stiffness and friction on the plastic 
behavior of a granular packing undergoing ratcheting will 
be also investigated in the following sections. 



V. HYSTERETICAL BEHAVIOR 

History dependence is one of the most essential fea- 
tures of granular soils. In our simple model, we have 
shown the existence of hysteresis both in the shakedown 
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FIG. 6: Evolution of the resilient parameters with the number 
of cycles: Resilient modulus Mr (top) and Poisson ratio £ 
(bottom) . The curves show the measures of these magnitudes 
for different values of the deviatoric stress Act. The data 
in the figure correspond to the simulation of a system with 
friction coefficient fx = 0.1, normal stiffness k n — 2 • 10 6 N/m, 
normal damping l/j n = 4 • 10 2 t s , and tangential damping 
1/7* = 8 • 10H B . The confining pressure is Pq — 6-10 k n . 



and in the ratcheting regime. This has forced us to iden- 
tify two different components to the total strain, namely 
the permanent and the resilient strain. In any stress cy- 
cle, the sliding contacts behave differently in the loading 
and un-loading phase, leading to a different stiffness of 
the material in each of these phases. In this section, we 
are interested in the shape of the cycles and, more specif- 
ically, in its relationship with the evolution of the area 
closed by the strain-stress loop. If we assumed that the 
deformation in both spatial directions is approximately 
the same, this area is the dissipated energy within the cy- 
cle. This energy relaxes during the post compaction from 
an initial high value to a constant value [23], reflecting 
the similarity of the hysteresis loops in the ratcheting 
regime (see Figure 1). This final value is plotted in Fig- 
ure 7 for different deviatoric stress. A clear power law 
behavior is observed in a wide range of values above the 
shakedown regime. 



For the purposes that will be seen next, let us introduce 
the following variables: 



7 = 7o + -y - 7, 
„ _ Act ct 2 - 01 
~ ~2 ~Pfa ' 



(16) 
(17) 



Being 70 the permanent strain accumulated up to the 
end of the previous cycle we are interested in. 

We express in Figure 8 the limit cycle on Figure 1 
on these new variables. The best-fit curve to the points 
in the loading and unloading are also included. These 
curves can be expressed, using the scaled variables: 



1 




in the loading. And 



1 * R ( A(J 

^=M~ R q - Bu ({- 



_*2 



(18) 



(19) 



in the unloading phase. Bl and Bjj are positive con- 
stants dependent on the confining pressure, but indepen- 
dent on the maximum deviatoric stress (Act). Note the 
use of the resilient parameter Mr in the previous expres- 
sions. From these formulas, it is then trivial to find the 
area of the cycle (Ah)- 



ih 



(B L + Bu) 



Act 



02-01 

Po 

Aer/2 



d-f 



(20) 
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Act 
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-Act j . 



Due to our definition of q* and 7*, the area Ah is in 
fact the same as the area enclosed by the stress-strain 
cycle in Figure 1. Our simple calculation explains why, 
given the nature of the stress-strain cycles obtained in 
our model, the power law behavior on Figure 7 should 
be expected. The explanation shown here somehow re- 
sembles the Rayleigh law for magnetization of ferromag- 
netic materials under low inductions [48]. Also in this 
case, the hysteresis energy loss (the area of the induction 
versus magnetization loop) behaves like the cube of the 
induction. This power law in ferromagnetic materials 
results from the quadratic dependence of the magnetic 
field on the magnetization. This is analogous to Eqs. 
( 19) and ( 18) except for the fact that Bl ^ Bu, which 
reflects the asymmetry of the loops in the granular ratch- 
eting regime. It is interesting to observe that the power 
law is identical to the one found for the dependence of 
the strain-rate on the deviatoric strain, as shown in the 
previous section. In fact, the closed-loop approximation 
given by Eqs. ( 19) and ( 18) is not strictly valid in 
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the limit q* — > 0. The error of this quadratic approxi- 
mation is of the order of 0(cr 3 ), and must be related to 
the cubic dependence of the strain accumulation on the 
load amplitude. A micro-mechanical explanation of this 
Rayleigh-like law in granular ratcheting is still an open 
issue. 



their dependencies on the model parameters will help to 
gain insight into the overall plastic response of the ma- 
terial. 



VI. PERMANENT STRAIN ACCUMULATION 
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FIG. 7: Variation of the area enclosed by the stress-strain 
cycle Ah, for different values of Act. The area is scaled with 
the confining pressure. The dashed line shows the power law 
y oc x 3 . The data in the figure correspond to the simulation 
of a system with friction coefficient fi — 0.1, normal stiff- 
ness k n = 1.6 • 10 6 N/m, tangential stiffness k t = 0.33fc„, and 
normal damping l/j n = 4 • 10 3 t s . The confining pressure is 
Po = 6 ■ !CT 3 fc„ and the damping coefficient j t = 8t s . 
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FIG. 8: Hysteresis stress-strain loop in the new variables 7* 
and q* . The solid points are the result of the simulation shown 
in Figure 1 (N — 1000). The solid lines are the best- fit to 
the expressions ( 19) and ( 18). The values of the constants 
for the theoretical lines are Bl = 0.04543, and Bjj — 0.05554. 



In the ratcheting regime the factors follow Bjj > Bl- 
It is still to be determined which precise effect has the 
behavior of the sliding contacts on this observation. A 
better understanding of the nature of these constants and 



The influence of macro-mechanical magnitudes and the 
microscopic parameters of the model on the accumula- 
tion of permanent strain will be shown in this section. 
This will be done by measuring the strain rate in simu- 
lations were the confining pressure, the deviatoric stress, 
the friction coefficient or the stiffness of the contacts are 
changed, while the rest of the parameters are kept fixed. 



A. Influence of the confining pressure and 
deviatoric stress 

Among all the possible parameters affecting the plas- 
tic behavior of a granular sample, the dependence on the 
confining pressure and on the deviatoric stress are known 
to be the most relevant ones [8]. Since P is measured 
in units of the normal stiffness, Po = Pok n , in our sim- 
ple model, there are two equivalent ways of studying the 
effect of the confining pressure: On the one hand, the nor- 
mal stiffness of the contact can be changed while main- 
taining the ratio kt/k n constant. On the other hand, the 
effective pressure Po can be increased. In order to inves- 
tigate the importance of the stress history of the sample, 
both methods have been used and the results are shown 
on Figure 9. (a). In each of the simulations, the system 
was first homogeneously compressed, and then subjected 
to cyclic loading. A power law relating the change of 
strain per cycle, A7/A./V, to Po/k n in a wide range of 
values is found in our simulations. The best fit of the 
points leads to the linear behavior: 



A7 P, 
AN " k n 



(21) 



Dispersion of the data with respect to the empirical law 
in Eq. ( 21), is a direct consequence of the dependence 
of the final strain rate on the preparation of the mate- 
rial. Different confining pressures imply a different post- 
compaction process [23] and therefore a different density 
of the sample before cyclic loading. The range of den- 
sities involved in Figure 9. (a) goes from solid fractions 
<F = 0.82 to $ = 0.9. Our results show, in fact, that the 
strain-rate seems to be much more sensitive to changes 
in the density than the resilient parameters. This makes 
the investigation of the strain accumulation more diffi- 
cult, limiting also the accuracy of our results on the rela- 
tionship between the basic parameters of the system and 
the strain-rate. 

This history dependence of the material is not observed 
in part (b) of Figure 9, where the strain-rate accumula- 
tion is plotted versus the deviatoric stress for the same 
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FIG. 9: Strain-rate dependence on the confining pressure, Po, 
and the deviatoric stress Act. The solid line represents the 
best-fit power law. The simulation details are those of Fig- 
ure 7. Data on the top graph correspond to Ac = 0.2 and 
tangential damping 1/74 = 8 ■ W 2 t s . Solid circles were ob- 
tained keeping k n constant and varying Po. The open circles, 
on the contrary, are the result of a series of simulations in 
which k n was changed. The solid line on this graph shows a 
linear behavior. Data for the plot on the bottom correspond 
to Po = 6 ■ 10~ 3 fc n and y t = 8t s . The solid line represents the 
power law y oc x 3 . This is close to the power law fitting in 
polygonal packing, whose exponent lies between 2.7 and 2.9 
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FIG. 10: Dependence of the strain rate on the friction coeffi- 
cient (i. The data in the figure correspond to the simulation 
of a system normal stiffness k n = 1.6 • 10 6 iV/m, tangential 
stiffness k t = 0.33k„, normal damping l/j n = 4 • 10 3 i s , and 
tangential damping 1/74 = 8t s . The stress conditions are 
Po = 10~ 3 • k n and Aa — 0.1. The solid fraction of the initial 
condition is $ = 0.93. The solid line shows the law y = x~ 2 . 



systematic ratcheting effect can be found. For the pa- 
rameters used in the simulation shown in Figure 10, this 
limit value is \i = 0.05. The strain-rate is maximal at this 
friction, and (as observed in the figure) the strain-rate 
decreases from this point, as friction is increased. The 
explicit dependence on the friction coefficient follows the 
power law: 



A 7 , V. 



2. ±0.05 



(22) 



Figure 11 shows the variation of the permanent strain 
accumulation rate with the stiffness ratio for different 
samples prepared with the same confining pressure Po 
and normal stiffness k n . A power law behavior with a 
negative exponent is found. The best fit of the points of 
the figure gives: 



initial configuration of disks with solid fraction <f> = 0.85. 
The measures indicate a clear potential dependence of the 
strain-rate with Ac. Also an exponential behavior (with 
exponent m = 2.8±0.1) has been reported in a polygonal 
packing [1]. 

B. Influence of the micro-mechanical parameters 

The strain-rate behavior as friction changes is slightly 
more complicated, if compared to the other parameters 
studied. For very low friction, no ratcheting is observed 
in the sample. Above a certain value of ji, however, a 



A7 (kt\ , x 

indicating that stronger tangential forces produce a 
higher rate of the deformation. 

An interpretation of these power law relation could be 
done by exploring the statistical distribution of the con- 
tact forces and its evolution during the loading stage. An 
important parameter is the mobilized angle a = \ft\/f n , 
which is bounded by the sliding condition a = /1. The 
statistical distribution of this variable is rather constant 
except for a peak at \x given by the sliding condition. The 
value of this peak depends on the friction coefficient. For 
small values of ji a large number of contacts can reach 
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FIG. 11: Dependence of the strain-rate on the stiffness ratio 
k t /k„. Data correspond to the simulation of a system with 
normal damping l/j n = 4 • 10 3 i s , and tangential damping 
1/jt = 8 ■ 10 2 i s , solid fraction <E> = 0.845 ± 0.005, and friction 
coefficient /i = 0.1. The stress conditions are kept constant, 
Po = 10~ • k n and Act = 0.2. The solid line represents the 
power law y oc a; -0 ' 3 . 

the sliding condition so that the ratcheting response is 
expected to be large. For big values of \i only a few num- 
ber of contacts can reach the sliding conditions, which 
produces a small ratcheting response. A quantitative 
explanation for the power law dependence will require 
to calculate the evolution of the statistics of the sliding 
contacts and the contribution of the sliding to the global 
dissipation, but this is beyond the scope of this work. 

VII. RESILIENT RESPONSE 

Most theoretical models for the resilient response are 
based on curve fitting procedures, using data from biaxial 
or triaxial tests. One of the most popular and earlier 
models is the so-called k — 9 model [35], in which the 
resilient modulus is supposed to depend only on the mean 
stress 9: 

M r (0) = k(^j , (24) 

where k and n are material constants, r] is a universal 
constant in units of stress (included for normalization), 
and 9 is the absolute value of the first invariant of the 
stress tensor: 

= |tr(*)l- (25) 

Many alternatives to and modifications of this model 
have been introduced, which are extensively used in prac- 
tice [34, 49, 50]. One of the main restrictions of the k — 9 
model is the assumption of a constant Poisson ratio. Sev- 
eral studies have shown that the Poisson ratio is not a 



constant in the granular case, but varies with the applied 
stresses [36] . Another drawback of the model is that the 
effect of the deviatoric stresses on the resilient modulus 
is neglected. A straightforward modification of the k — 9 
model accounting for this latter restriction reads [37]: 

M r (0,Aa) = k(^y '{^Y '. (26) 

Note that, with respect to equation (24), a new material 
constant m has been introduced. In the simplest approx- 
imation both exponents are assumed identical n = m 
[38]. 

The validity of the k — 9 model will be checked in this 
section. Note that, in the case of cyclic loading, given a 
fixed Act, the dependence of the resilient modulus on 9 is 
similar to its dependence on P . Results will be shown on 
the influence of the confining stress and deviatoric stress 
on the resilient modulus and Poisson ratio. In the latter 
case, it will be particularly interesting to investigate the 
limit of validity of the common assumption of a constant 
Poisson ratio for granular matter. 

A. Influence of the confining pressure 

Figure 12 indicates that the k—9 model is in fact a very 
good approximation in the ratcheting regime for a wide 
range of pressures of Po. The best fit to the empirical 
law of Eq. ( 24), gives n = 0.34 ±0.02. This value agrees 
well with the experimental values in [36], where results 
on gravel show a power law with exponent n = 0.31. 

The Poisson ratio behaves in a completely different 
way. For low pressures, it decreases gradually as the 
pressure becomes higher. For P > 0.01fc„, however, 
there is a change on the trend, and ( grows fast with 
Po. This reflects a higher anisotropy of the deviatoric 
strain in systems compressed under a high pressure. Nev- 
ertheless, our results justify the use of a constant value 
of ( in a first approximation, for a wide range of P , 
10~ 4: k n < Pq < 10~ 2 fc„. The most common estimate 
(C = 0.35), however, slightly overestimates the values 
obtained in most of our simulations. 

B. Influence of the deviatoric stress 

Two stages are clearly distinguished in the behavior 
of the resilient parameters as a function of Act. For low 
values of the deviatoric stress, close to the shakedown 
regime, the resilient parameters remain approximately 
constant. Poisson ratio, remains closer to the indicated 
value C ~ 0.35 which is the empirical fixed value usually 
assumed for unbound granular matter [36]. This value 
is shown in Figure 13 with a solid line. For Act > 0.1, 
however, £ shows a strong dependence on the deviatoric 
stress Act. 

A simple empirical polynomial law is proposed in ref- 
erence [36] for the dependence of ( on the ratio of the 



11 




0.001 



le-04 le-03 le-02 le-01 




0.01 



le-05 le-04 



le-03 le-02 
P /k n 



le-01 le+00 



FIG. 12: Variation of the resilient parameters with the con- 
fining pressure Pq: resilient modulus Mr (top) and Poisson 
ratio ( (bottom). The conditions of the simulation are the 
same as in Fig. 9. The line in the left plot is the best fit to 
the k — 9 model. The solid line in the right figure is the value 
£ = 0.35, estimation for the Poisson ratio of granular mate- 
rials. The different symbols refer to two different methods 
explained in the text to study the influence of the confining 
pressure on the system. 
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FIG. 13: Variation of the resilient parameters with the loading 
intensity Act. The simulation details are similar to those in 
Figure 9 but with k n = 2 • 10 6 iV/m, and P = 6 ■ 10~ 4 k n . 
The best fit curve to a second order polynomial is plotted for 
the values of Mr in the top graph. In the bottom (Poisson 
ratio), the solid line corresponds to the value £ = 0.35 and 
the dotted line to the best fit to equation y(x) = a + bx + cx 2 
(details are given in the text). 



deviatoric and volumetric stresses. Although the range 
of values studied in this experiment is larger than the 
one presented here, our results confirm that the values of 
the Poisson ratio follow a second order polynomial law 
on Act, being the best-fit curve C = 0.336(±0.001) - 
0.208(±0.001)Acr + 3.061(±0.001)(Acr) 2 . This curve is 
plotted in the lower part of Figure 13. 

As opposed to the behavior of Poisson ratio, the re- 
silient modulus decreases as Act increases, the depen- 
dence is also polynomial. In Figure 13 (top), the curve 
y(x) = 335.7-316.8a; + 229.1x 2 is plotted. Note that this 
result disagrees with the simplification of the generalized 
k — model (m = n) of equation (26). The general law 
seems to be a better approximation in a wide range of 
values of the deviatoric stress, where the system shows 
neither collapse nor shakedown. 

The dependence of the resilient parameters on the de- 



viatoric stress results from the anisotropy induced in the 
contact network for large deviatoric loads. Near failure, 
a significant number of contacts are open in the perpen- 
dicular direction of the load, resulting in a decrease of 
the stiffness as shown in the top of Figure 13. The in- 
crease of the Poisson ratio in the bottom of this figure 
is consequence of the formation of force chains, which 
enhance the anisotropy and leads to an increase of the 
effective Poisson ratio. A detailed description of the ef- 
fect of these force chains in the resilient response would 
require a detailed evaluation of the relation between the 
anisotropy of the contact network and the parameters of 
the anisotropic elasticity via fabric tensors [17, 51]. 
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FIG. 14: Variation of the resilient modulus with the static 
friction coefficient /i. The conditions of the simulation are 
the same as in Fig. 10. 



C. Influence of the micro- mechanical parameters 

Figure 14 shows the change of the resilient modulus 
with friction. Mr grows for small frictions. However, 
the curve seems to reach a saturation level for frictions 
H « 0.4. 

Changing the ratio of contact stiffness (Fig. 15), a 
power law dependence of Mr is observed for k n /k t < 0.1, 

fu \ - 28 

M R oc f fM , being the exponent 0.28 ± 0.03. For 

stiffness ratios closer to unity k t /k n w 1, the resilient 
modulus remains approximately constant or even de- 
creases. The Poisson ratio also appears to be constant 
for k t < 10~ 3 • k n . Above k t /k n = 0.001, ( decreases to 
values below the reference value C = 0.35. For k t > k n , 
C starts growing again. 
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FIG. 15: Influence of the ratio of contact stiffness k t /k n on 
the resilient parameters. The details of the simulation are 
those of Fig. 11. The solid line shows a power law with 
exponent 0.28 in the top. The one at the bottom marks the 
value £ = 0.35. 



VIII. DISCUSSION AND FINAL REMARKS 

A characterization of the material response in the gran- 
ular ratcheting has been presented in terms of the strain- 
rate, resilient modulus and the Poisson ratio. Studying 
the dependence of these parameters on the conditions of 
the biaxial test (stress configuration) and the main mi- 
croscopical constants of the sample (friction and contact 
stiffness) we confirmed the persistence of the granular 
ratcheting in many different conditions and systems. 

Given a compressed sample subjected to a biaxial test 
in which a cyclic loading is switched on, the system 
adapts to the new situation accumulating deformation 
and dissipating energy at a relatively high rate. After 
this post- compaction stage, the dissipated energy, both 
resilient moduli and the strain-rate reach stationary val- 
ues. The duration of the adaptation stage basically de- 
pends on the deviatoric stress, and is usually shorter for 
the resilient moduli than for the strain-rate [23]. If the 



deviatoric stress is small enough, the perturbation intro- 
duced by the cyclic loading shakes down. The material 
adapts to the new situation so that there is no further 
accumulation of permanent strain. Above this limit the 
material accumulates a certain amount of strain in each 
cycle. If the stress is below the collapse limit, the per- 
manent strain accumulated after each cycle is constant. 
This is the so-called granular ratcheting, which has been 
described both experimentally [15, 52] and in simulations 
[1, 23]. 

Identical repetition of the strain-stress cycles is among 
the main characteristics of the granular ratcheting. This 
periodicity reflects the weak dependency of the resilient 
moduli on the stress history and, in the particular case 
of cyclic loading, on the number of applied cycles [8] . In 
all the simulations, a steady and stable resilient response 
is reached after some initial cycles. This kind of sim- 
ple behavior is expected as long as the applied deviatoric 
stress remains below the collapse limit. Although many 
factors may influence the plastic response of the system, 
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there is a simple characterization of the deformation in 
the ratcheting regime, in terms of the strain-rate and the 
resilient moduli. This description takes advantage of the 
empirical fact that these magnitudes do not change in 
the ratcheting regime. We have investigated both micro- 
mechanical and macro-mechanical factors influencing the 
plastic response of the material, i.e. the dependency on 
the number of cycles, static friction, the confining pres- 
sure, the deviatoric stress and the stiffness. 

It was shown that the use of a constant Poisson ratio 
is a good approximation in most cases. It seems to be 
unsuitable, however, for very high confining pressures, 
very high deviatoric stresses, or for low values of the fric- 
tion coefficient. The value for £ estimated through our 
simulations would be slightly below the empirical value 
0.35, assumed in many models of the resilient response 
of granular materials. This might be a consequence of 
the simplicity of the visco-elastic model, which does not 
include all the mechanisms involved in a real biaxial ex- 
periment. 

M R is a measure of the macroscopical stiffness of the 
material. Our results show that it is higher for strongly 
frictional materials. We also found that although prepar- 
ing the sample with a higher confining pressure increases 
its stiffness, increasing the deviatoric stress reduces the 
stiffness of the packing. 

Both the strain-rate and the resilient modulus Mr 
show a power law dependence with the confining pres- 
sure and the ratio of contact stiffness. The power law 
is similar for both magnitudes in the case of the confin- 
ing pressure, but they have an opposed dependence on 
k t /k n . The dependence of Mr on the deviatoric stress 
is a second order polynomial. The generalization of the 
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